#include<iostream>
#include<stdio.h>
#include<vector>
#include<math.h>
#include<string.h>

using namespace std;
typedef double Type;
const int N=2000;

void Cholesky(Type A[][N],Type L[][N],const vector<Type> B,int n);
vector<Type> Solve(Type L[][N],vector<Type> X,int n);

void Cholesky(Type A[][N],Type L[][N],const vector<Type> B,int n)
{
  for(int k=0;k<n;k++)
    {
      Type sum=0;
      for(int i=0;i<k;i++)
	sum += L[k][i]*L[k][i];
      sum=A[k][k]-sum;
      L[k][k]=sqrt(sum>0?sum:0);
      for(int i=k+1;i<n;i++)
	{
	  sum=0;
	  for(int j=0;j<k;j++)
	    sum += L[i][j]*L[k][j];
	  L[i][k]=(A[i][k]-sum)/L[k][k];
	}
      for(int j=0;j<k;j++)
	L[j][k]=0;
    }
  /*output L*/ 
  for(int i=0;i<n;i++)
    {
      for(int j=0;j<n;j++)
	cout<<L[i][j]<<" ";
      cout<<endl;
    }
  /*output the result*/
  cout<<endl;
  vector<Type> X=Solve(L,B,n);
  vector<Type>::iterator it;
  for(it=X.begin();it != X.end();it++)
    cout<<*it<<" ";
  cout<<endl;
}

vector<Type> Solve(Type L[][N],vector<Type> X,int n)
{
  /*LY=B => Y*/
  for(int k=0;k<n;k++)
    {
      for(int i=0;i<k;i++)
	X[k]-=X[i]*L[k][i];
      X[k]/=L[k][k];
    }
  /*L^tX=Y => X*/
  for(int k=n-1;k>=0;k--)
    {
      for(int i=k+1;i<n;i++)
	X[k]-=X[i]*L[i][k];
      X[k]/=L[k][k];
    }
  return X;
}
